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Abstract. 

Realistic, first-principles-based interatomic potentials have been used in 
molecular dynamics simulations to study the effect of cation composition 
on the ionic conductivity in the Zr2Y207-Y3Nb07 system and to link 
the dynamical properties to the degree of lattice disorder. Across 
the composition range, this system retains a disordered fluorite crystal 
structure and the vacancy concentration is constant. The observed trends 
of decreasing conductivity and increasing disorder with increasing Nb^"*" 
content were reproduced in simulations with the cations randomly assigned 
to positions on the cation sublattico. The trends wore traced to the 
influences of the cation charges and relative sizes and their effect on 
vacancy ordering by carrying out additional calculations in which, for 
example, the charges of the cations were equalised. The simulations 
did not, however, reproduce all the observed properties, particularly for 
YaNbOr. Its conductivity was significantly overestimated and prominent 
diffuse scattering features observed in small area electron diffraction studies 
were not always reproduced. Consideration of these deficiencies led to a 
preliminary attempt to characterise the consctjuencc of partially ordering 
the cations on their lattice, which significantly affects the propensity for 
vacancy ordering. The extent and consequences of cation ordering seem to 
be much less pronounced on the Zr2Y207 side of the composition range. 
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1. Introduction 

Oxide-ion conducting materials are the subject of intensive research activity due to 
their technological apphcations in Sohd Oxide Fuel Cells (SOFCs), oxygen separation 
membranes and gas sensors. Cubic, fluorite-structured, binary compounds of stoi- 
chiometry AO2 are of particular interest, especially when some of the host cations 
are replaced by species of a lower valence to produce anion-deficient phases like 
A'^'^^B^~^02-x/2- These compounds have charge compensating vacancies which al- 
low the remaining anions to move more rapidly throughout the crystal and this leads 
to impressive values of the ionic conductivity (cr > 0.1 il^^cm^^, at T = 1500 K) as 
observed in Af203-doped Zr02 (M = Sc or Y) and M203-doped Ce02 {M = Y or Gd). 

At the simplest level the ionic conductivity would be expected to increase with 
the number of vacancies, but this is not the case. In Zr^^^Y^"'"02_a/2 the conduc- 
tivity decreases with increasing x for x > 0.16 and, furthermore, for a given x the 
value of the conductivity depends on the identity of the dopant cation [T] . For exam- 
ple, Zr^t^Sci^"'"02-a;/2 is a better conductor than Zr^ta:YiJ'''02-a;/2 These effects 
are due to a reduced mobility associated with the ordering of the vacancies, which 
may be influenced by vacancy-vacancy, vacancy-cation and cation-cation interactions. 
Bogicevic and co-workers [3] examined these effects in stabilized zirconias by using 
first-principles, zero temperature electronic structure calculations. They found that 
vacancy ordering is most strongly governed by vacancy-vacancy interactions, followed 
by vacancy-dopant and, weakest of all, dopant-dopant interactions. The influence of 
the identity of the dopant could be understood in terms of a balance between com- 
peting electrostatic and elastic effects, which are associated with differences in the 
cation charges and radii. They showed that the best way to minimize vacancy-dopant 
association is to choose dopants where the strain term (i.e. the dopant ionic radius 
relative to that of the host cation) counteracts the electrostatic interaction, rather 
than simply matching the ionic radii. Some of these factors have also been identified 
in other studies [H [5l [6l [7| . 

We have carried out Molecular Dynamics (MD) simulation studies of a similar 
class of materials with realistic first-principles based interaction models in order to 
examine the influence of these effects at the experimentally relevant temperatures (> 
1000 K). The interaction potentials are capable of accurately reproducing experimental 
conductivities etc., and long range effects of strain are incorporated by our use of 
considerably larger simulation cells (typically 4x4x4 unit cells but up to 10 x 
10 X 10) than are accessible in direct first-principles calculations. In the present 
work we focus on the ternary system Zro.5_o.52:Yo.5+o.252:Nbo.25a;Oi.75 which retains 
a disordered fluorite structure with a constant vacancy count as x is varied (one in 
8 oxide ions is "missing" ) , which allows us to selectively study the effect of varying 
dopant cation charge and ionic radius (Nb^"*" is significantly smaller than Y'^+ since 
R^^^ = 0.74, 0.84, 1.019 A for , Zr*+ and Y'''+ respectively 0). Although 

the number of vacancies is fixed across the whole series, the conductivity changes by 
almost two orders of magnitude from a; = to x = 1 at 1000 K [10 . Furthermore, 
it was observed that the intensity of diffuse scattering, related to the lattice disorder, 
decreases as the conductivity increases which presented the prospect of linking 
structural (crystallographic) studies to the effects on the conduction mechanism and 
the conductivity. 
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Consequently, we have examined [HI [11] this system with a combination of 
impedance spectroscopy, powder neutron diffraction - including the analysis of the 
total scattering using bond valence sum constrained Reverse Monte Carlo modelling 
(RMC) [H [HI HH - and MD simulations. The RMC analysis indicated that real- 
space configurations based on a fluorite structure with the dopant cations randomly 
arranged within the cation sublattice, i.e. a disordered-fluorite (rf-fluorite) structure, 
gave an excellent representation of the diffraction data across the composition range. 
The MD simulations were initiated from such d-fluorite configurations. They accu- 
rately reproduced the average structure extracted from the diffraction data using RMC 
at the level of the partial radial distribution and bond-angle distribution functions. 
Furthermore, the simulations reproduced the temperature-dependence of the experi- 
mental conductivity for Zr2Y207 {i.e. the x—Q system) and the trend of decreasing 
conductivity with increasing Nb content across the series. A substantial part of the 
present paper will be concerned with a detailed analysis of these simulations and exam- 
ination of the chemical factors which contribute to this trend. These analyses enable 
us to track the competition between Coulomb and strain effects on the vacancy or- 
dering in the d-fluorite structure and to explain the consequences for the ionic mobility. 

Other considerations suggest, however, that these simulations may miss an impor- 
tant contributor to the conductivity of these materials. Small area electron diffraction 
(SAED) studies have been carried out on these (and related) materials by several 
groups [m [16j [17] . This technique gives single-crystal diffraction information. The 
SAED studies show sharp diffuse scattering features at particular positions in recip- 
rocal space which are the result of a structural modulation which is not intrinsic to 
the d-fluorite structure [T5J [THl [13 • The diffuse peaks are associated with a pattern 
of strain caused by vacancy ordering and, according to some authors [TS] , to partial 
ordering of the dopant cations. The ordering is of intermediate range, one analysis 
suggesting a correlation length of ~22 A in YaNbOy [15]. The diffuse peaks suggest 
a structure with a local ordering tendency in which the anion vacancies are aligned 
in pairs along the (111) directions within the x —1 (YsNbOy) composition, in a man- 
ner related to the pyrochlore structure. As a: — > (Zr2Y207) the diffuse scattering 
changes and, as interpreted by Irvine et al |16j . is associated with vacancy alignment 
along (110) directions, which resembles the situation within the C-type structure of 
Y2O3. According to these authors, the decrease in conductivity as x -^1 is caused 
by the effects of these different patterns of vacancy ordering on the anion mobility, 
though no clear explanation of this structure-property relationship has been provided 
to date. To anticipate the story below somewhat, not all of these peaks are seen in 
diffraction patterns calculated from our d-fluorite simulations, even if they are initi- 
ated with configurations in which the oxygen ions are placed in positions consistent 
with the pyrochlore vacancy ordering. We have, therefore been led to perform a new 
set of simulations in which we explore the effect of some degree of local cation ordering. 

In this paper, we begin with a detailed examination of the structure and 
ion dynamics in the d-fluorite simulations which successfully reproduce the trends 
observed in these properties in the Zro.5_o.5a;Yo.5+o.25a:Nbo.25a:Oi.75 series. In order 
to clarify how they are affected by the cation properties, we consider the effects of 
modifying our first-principles interaction potentials by equalising the cation charges 
and the range of the short-range repulsive interactions between the cations and the 
oxide ions. We then consider the structure of the material from the perspective of the 
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vacancies and demonstrate similar cation-specific vacancy-ordering tendencies to those 
which have been seen in the ab initio studies of stabihsed zirconias [3J . For the reasons 
described above, we are then led to examine the consequences of partially ordering 
the cations for the energetics, the diffuse scattering and the conductivity. The studies 
here are less definitive than the d-fluorite ones, as we cannot determine the degree 
of ordering within our MD simulations and must rely on drawing conclusions from 
different postulates for the order. The partially ordered structures necessarily have 
a lower entropy than the fully disordered ones, and therefore become unstable with 
respect to them at high temperatures. However it is only at very high temperatures 
(where these phases are unstable) that the cations move sufficiently to allow the 
structure to equilibrate. 

2. Disorder and mobility in the rf-fluorite structure of Zr2Y2 07 and 
YgNbOr 

The simulation methodology is summarised in [Appendix A[ The simulations are thor- 
oughly equilibrated at high temperatures (1500-2000 K) in constant pressure sim- 
ulations before obtaining dynamical information from constant volume runs at the 
zero- pressure volume. At these high temperatures the oxide ions are diffusing on the 
simulation timescale but there is no significant exchange of cations. We also quenched 
down to room temperature and conducted further short runs in order to examine 
structural information and compare with diffraction data collected at the same tem- 
perature. Note that, because the ionic mobility drops sharply as the temperature is 
reduced during the quench, that these runs are not fully equilibrated. To some extent, 
this will also be true of experimental samples which are prepared at high temperature 
and quenched, though the rate of cooling is much slower than for the simulated sys- 
tems. 

We begin by discussing results on the d-fluorite structure, which are initialised 
by randomly distributing the cation species over the cation sublattice and randomly 
assigning oxide anions to the anion lattice. Theoretical calculations show the energy 
scale associated with cation ordering is very small O [TS] in yttria-stabilised zirconia 
and, as we showed in the previous paper [8], extensive powder diffraction and RMC 
analyses across the whole series Zro.5_o.5a;Yo.5+o.25a;Nbo.25a;Oi.75 provided no evidence 
for cation ordering, in agreement with previous studies [10]. We also demonstrated 
the excellent agreement between the present simulations and the RMC data, as well 
as a very good reproduction of the temperature dependence of the experimental con- 
ductivity for Zr2Y207. In the present paper we will only discuss results for the end 
members (a;=0 and x=\) of the series, i.e. 2,^2^207 and YaNbOy. 



2.1. Cation-induced disorder 

The disordered character of the anion sublattice can be appreciated from the oxygen- 
oxygen radial distribution functions taken from the room temperature runs, .go-o('') 
(see figure [TJ. They are compared with those expected from a perfect fluorite structure 
with the same lattice parameter (and with the 5o-o('') obtained from the RMC anal- 
ysis of the experimental data [8J). The first peak in go-o{r) for YaNbOy is shifted by 
0.13 A from the expected position for an ideal fluorite structure and the second peak 
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(which indicates the distance between two next-nearest neighboring oxygen ions) is 
broadened and presents a number of subsidiary features. Zr2Y207 shows similar but 
less marked deviations from the ideal fluorite order, as already noted from the diffuse 
scattering intensity P HU], 




Distance / A 



Figure 1. Comparison between the simulated (solid line) go-o(r) (obtained 
from the QUAIM potential described m ^81) in Zr2Y207 (bottom) and 
YsNbOj (top) and the one obtained from a perfect fluonte structure (dashed 
Itne) at room temperature. The experimental data (dots) from [SJ obtained 
via bond valance sum constrained RMC are shown as well. 

In the fluorite structure, the anions sit in a tetrahedral site. Because the cation 
species are distributed randomly in the d-fluorite Zr2Y207-Y3Nb07 system, the ver- 
tices of these tetrahedra may be occupied by different combinations of cations. A 
snapshot of the simulation cell in YaNbOy can be seen in figure [2] The two black 
circles highlight those oxide ions that are surrounded by Y^+ cations only and the 
red squares those with a significant number of Nb^"*" ions in their coordination en- 
vironment. Notice that the former are much more tightly grouped about the lattice 
site than the latter. The notion that the Nb^"*" cations are somehow responsible for 
the disorder observed in YaNbOy is confirmed by the fact that the Y-0 and Nb-0 
distances, obtained from the RMC analysis of the neutron scattering data in reference 
[S] (or analogously from our simulations), are quite different (dNb-o — 1-96 A dY-o= 
2.28 A). In a perfect fluorite structure, the anion-cation distance is dan-cat — \/3/4 
a = 2.27 A, where a is the lattice parameter obtained from the RMC analysis. It is 
therefore evident that the Nb"*"^ ions are attracting the oxygen ions very tightly, thus 
disrupting the perfect fluorite order. This effect is less signiflcant in Zr2Y207 because 
the cation charges are 3-|- and 4+ in this case and because the cation radii are more 
similar. 
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Figure 2. Snapshot along ttie z axis of the simulation box m Y^NbOj. The 
color code is whtte, green and red for Yttrium, Niobtum and Oxygen ions 
respectively. (In black and white, the colour code is white, light grey and 
dark grey for Yttrium, Niobium and Oxygen ions respectively.) The two 
black circles show those oxygen ions which are surrounded by Yttrium ions 
only, whilst the two red squares show those which are surrounded by both 
Niobium and Yttrium ions. 



2.2. Cation effects on anion mobility 

The inhomogeneities induced by this strain lead to a significant difference in conduc- 
tion mechanisms between Zr2Y207 and YsNbOy. In figure |3]we show the probabiHty 
distribution (calculated as in [8 ) for the individual oxide ions to make a certain num- 
ber of diffusive hops (by the 0-0 nearest-neighbour separation) during the course of 
the simulations at 1500 K. In YaNbOy a large fraction of anions do not hop at all 
during the whole simulation whereas a few anions possess a high mobility. The for- 
mer make no contribution to the conductivity. Zr2Y207, on the other hand, has a 
more homogeneous behaviour with fewer immobile oxygen ions and a probability dis- 
tribution which peaks about the average mobility of 2-3, with relatively few outliers. 
Further examination (see below) shows that the immobile oxide ions in YaNbOr are 
in tetrahedral sites where all the vertices are occupied by Y'^+ cations. 



As a consequence, the conductivity in YaNbOr exhibits an "aging" effect, espe- 
cially at low temperatures. The average mobility shortly after initialising the simula- 
tion with randomly placed oxide ions is high, but drops over a substantial period of 
time (~1 ns) presumably because, on this time scale, oxide ions find initially empty 
Y3+_j-ich sites and become trapped there. Another consequence is that different ini- 
tial configurations give different mean diffusion coefficients even after very long runs. 
This is probably caused by the variation in the numbers of different types of tetra- 
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Figure 3. Number of ions with certain mobility against the mobility itself for 
YsNbOj (solid line), ZriY-i,Oi (dashed line) , equal-charge Y^NbOr (dotted 
line) and equal-radius Y^NbOr (dot-dashed line). 



hedral sites in the initialisations. Both of these effects make a rehable determination 
of the conductivity difficult for YaNbOy in the temperature range (< 1500 K) where 
experimental measurements have been made. Our best efforts give a conductivity 
(cAfB = 0.016 n~^cm~^ at 1500 K, where the subscript NE denotes Nernst-Einstein, 
see below) of about an order of magnitude larger than measured [Sj jlO,, . This value is 
much lower than the one obtained for Zr2Y207 {ctne — 0.08 ^l~^cm~^), which agrees 
with the experimental trend [S]- For Zr2Y207 neither of these effects is observed 
and reliable conductivity values can be obtained down to the highest temperatures at 
which measurements are available [|] In this case the slope of conductivity vs. inverse 
temperature agrees very well with experiment and the actual value at 1500 K is only 
slightly larger than experiment (by a factor of less than 2 for the data in |10J). 

The above-mentioned conductivity values were calculated from the ionic diffusion 
coefficients by using the Nernst-Einstein formula 

where c is the charge, p the number density and D the diffusion coefficient of the 
mobile species, ks is the Boltzmann constant and T is the temperature. The diffusion 
coefficients themselves have been calculated from the long-time slopes of a plot of the 
mean-squared displacement (msd) of individual ions versus time. The Nernst-Einstein 
formula assumes that the ions move independently, which can only be approximately 
true in reality. A better value for the conductivity is obtained, in principle, from an 

f There is a lower limit to the values of the diffusion coefficient and conductivity which may be 
determined from the mean-squared displacement of the ions as a function of time. To be reliable, the 
ions must on average be diffusing by more than a lattice parameter, and when the mobility is very 
low this requires the mean-square displacement to be measured over a very long time and eventually 
longer than we can afford to simulate. 
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integral over the charge current correlation function J(t)|19j: 

though this quantity has very poor statistics, especially when diffusion is slow. We 
were able to get a reliable value of the true conductivity from very long runs on the 
most highly conducting Zr2Y207 system. At 1500 K a = 0.05 il^^cm^^ which is 0.6 
times smaller than the Nernst-Einstein value and in much better agreement with the 
experimental values in [51 [TD] . At higher temperatures the differences between the two 
methods of calculating the conductivity become smaller. 



3. The effects of cation charges and sizes on the material's properties 

The simulations on the Zro.5_o.5a;Yo.5+o.252;Nbo.252;Oi.75 series in the rf-fluorite 
structure give the same trend of decreasing conductivity and increasing disorder with 
increasing x that is seen experimentally. In this section we will attempt to distinguish 
the effects of cation size and charge by comparing the results obtained with modified 
interaction potentials with the realistic models. 



3.1. The effect of the cation charge 

In order to examine the role of the high charge of the Nb^"*" cation in causing these 
effects we have carried out simulations on the x=l system in which both the "Nb" 
and "Y" ions are assigned a charge of 3.5+ (which maintains charge neutrality for 
this stoichiometry). We will call this system equal-charge YsNbOy. All other aspects 



of the interaction potentials were kept the same (see Appendix A I as in the runs de- 
scribed above. We started from a previous high temperature simulation on Y3Nb07 
and equalized the cation charges. All the simulations were equilibrated at constant 
pressure for 100,000 steps, during which the velocities were rescaled several times to 
keep the system at the required (high) temperatures for conductivity studies, and then 
longer (constant volume) runs were made to study the properties of this material. We 
also ran some shorter simulations at room temperature to study the static structure 
and obtain the radial distribution functions. 



In figure |4] we show the oxygen-oxygen radial distribution functions for room tem- 
perature equal-charge YsNbOy , YsNbOr itself, and a perfect fluorite structure with 
a matched lattice constant. It is clear that equalizing the cation charges increases the 
order of the oxide ion sublattice [f] The 5o-o('') for the equal-charge system shows 
a first peak close to the ideal fluorite position and a single prominent second peak, 
though still not as intense as in the ideal case. We also find that the cation-anion 
distances change significantly, dy-o — 2.23 A and dNb-o — 2.19 A. This is another 
sign that the system has become more ordered, as these become closer to the ideal 
cation-anion distance in this system, given, in this case, by \/3/4 ciequai-charge = 2.27 
A. A comparison between figure [Z] and figure [l] shows that equal-charge YsNbOy is 
also more ordered than Zr2Y207 which is consistent with the fact that Zr2Y207 still 
has two cations with different charges (4+ and 3-I-). 



§ A movie can be found in the online version of this paper which shows how YaNbOy becomes more 
ordered once the cation charges are equalized. 
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Figure 4. Comparison between the oxygen-oxygen g(r) for Y^NbO^ (red 
dotted line), equal-charge Y^NbOj (black solid line) and a perfect fluorite 
structure (red dashed line). 



Table 1. Diffusion coefficients. 



Material 


10" X Do (cm 


2 


YaNbOr 


0.924 




Equal-charge YsNbOy 


1.87 


T = 2000 K 


Zr2Y207 


2.87 




YgNbOr 


0.072 




Equal-charge YsNbO? 


0.21 




ZraYaOr 


0.44 


T = 1500 K 


Equal-radius YsNbO? 


0.37 




Equal-chargc&radius YaNbOr 


2.6 





Figure |3] shows that equahzing the cation charges (dotted hne) greatly reduces 
the number of immobile oxide ions, relative to Y3Nb07. In fact, the distribution 
comes to resemble that for Zr2Y207 quite closely and is indicative of a homogeneous 
pattern of diffusion. The mean squared displacement for the oxide ions in equal-charge 
YaNbOy was calculated at 1500 K and 2000 K and the respective diffusion coefficients 
are reported in table [T] From this it can be seen that equal-charge YaNbOy is more 
conducting than YsNbOy. Interestingly equal-charge YaNbOy has a lower diffusion 
coefficient than Zr2Y207, whereas we have just seen that the latter system has a more 
disordered oxide ion lattice. This suggests that there is some other factor at work 
besides the charge-induced strain in the lattice. 

In order to see how the dopant cation charge induces the disorder and affects the 
ionic mobility we calculated the mobility distribution as done in [5] but this time we 
also monitored the identities of the four cations surrounding the anion. As already 
anticipated above, the anions with no mobility (about 1/4 of the total number, see 
figure|3]) are in sites with a much higher Y''+ content than expected (« 3.65 as opposed 
to the average of 3 Y^+ and 1 Nb^+ which is expected from the YaNbOj stoichiometry) 
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Figure 5. Anion mobility as a function of the average tetratiedron 
volume for equal-charge (dashed line), equal-radius (dot-dashed line), equal- 
charge&radius Y^NbOr ( circles -i-line), Y^NbOj (solid hne) and Zr2Y2 07 
(dots). 

and the higher the mobihty the higher the content of Nb^+ ions in the tetrahedron 
surrounding the anion. This effect seems to be associated with the volume of the 
coordination tetrahedron (obtained from the volume of the tetrahedron with the four 
closest cations at its vertices). The Coulombic repulsion between two Nb^"*" ions is 
almost three times stronger than that between the Y'^^ cations and from the cation 
cation radial distribution functions we find that the nearest-neighbour cation-cation 
separations are dy-y = 3.69 A and (iNb-Nb = 3.86 A (the ideal fluorite value being 
a/-y/(2) = 3.78 A) so that the coordination tetrahedra are strained to an extent which 
depends on their Nb^"*" content with the Nb^"'"-rich ones significantly larger than the 
average. The anion mobility is plotted against the average volume of the surrounding 
tetrahedra in figure [s] (solid line). Alternatively, we could have plotted the mobility 
against the average number of Nb^"*" ions in the coordination shell and obtained a 
very similar curve. Note, firstly, that the tetrahedral volumes for YsNbOy span an 
appreciable range, from 6.4 to 7.5 A'^ and secondly, that the anion mobility is a strongly 
increasing function of the tetrahedron volume. This observation therefore explains 
the paradox that although highly charged Nb^"*" cations bind oxide ions strongly they 
promote local mobility but reduce the overall conductivity by straining the cation 
sub-lattice resulting in the formation of oxide ion trapping sites. We also carried out 
the same analysis for Zr2Y207 (dotted curve) and for equal-charge YaNbOy (dashed 
curve) and present the results in figure [s] In these cases, the range of tetrahedron 
volumes is much narrower than for YgNbOr (< 0.2 A^) suggesting a much less strained 
cation lattice and the dependence of the anion mobility on the tetrahedral volume is 
barely discernable. 

3.2. The effect of the cation size 

Proceeding by analogy with the equal- charge study, we can examine the effect of equal- 
izing the cation radii in YaNbOy. The short-range repulsion terms in our potentials 
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model the repulsion between two ions and therefore are intrinsically related to the ex- 
tent of the valence electron density around the ion. In figure |6] we show the short-range 
repulsion terms for the Y-0, Nb-0 and Zr-0 interactions (recall that these were gen- 
erated with the ah initio force-fitting strategy) . It can be seen that the range of these 
potentials is consistent with the fact that Y'^^ is the largest ion, Nb^"*" the smallest, 
with Zr^+ of intermediate size. To generate an equal-radius model for YsNbOy we 
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Figure 6. Short-range repulsion parts of the potential for the Y-0 (solid 
line), Nb-0 (dotted line) and Zr-0 (dot-dashed line) terms. 

used the Zr-0 short-range interaction parameters for the Y-0 and Nb-0 terms, but 
the Y'^'^ and Nb^"*" ions keep their formal charges ^ Simulations are run as previously 
described for the equal-charge case. 

The partial radial distribution functions for equal-radius YsNbOy looks quite sim- 
ilar to the one obtained for YsNbOr itself (figure [l] and reference [8]) and we have 
not shown them separately. Despite the equal radius potential construction, even 
the anion-cation distances are comparable (dNb-o = 1-93 A, dy-o = 2.29 A) to 
those of the parent compound. This confirms that the observed disorder is caused 
by the different cation charges in YaNbOy and that the different ranges of the repul- 
sive cation-anion interactions do not play an important role. Given this similarity in 
structure, it is not surprising that the equal-radius system shows a wide range of tetra- 
hedral volumes like Y3Nb07 itself and figure [s] (dot-dashed line) shows that, again, 
the oxide ion mobility tracks the tetrahedral volume like the parent compound. This 
leads to the inhomogeneous distribution of hopping probabilities, as shown in figure [3] 
(dot-dashed line) with a substantial fraction of immobile ions. 

In Table [l] we report the average diffusion coefficient obtained from a simulation 
at 1500 K on equal-radius YsNbOy and compare this value with equal-charge YaNbOy, 

II Note that the Y-O and Nb-O nearest-neighbour separations which result from these simulations will 
differ, because of the stronger Coulombic attraction to the Nb^+ ion. From an empirical perspective, 
i.e. the way that the ionic radius is normally assigned, the ions do not appear as of equal size; rather 
they have been treated as if the electron density distributions around the ions are of equal size in the 
potential model. 



Cation composition effects on oxide conductivity in the Zr2Y20-j-Y-^NhO'j system. 12 

YsNbOy and Zr2Y207 at the same temperature. Interestingly, the equal-radius ma- 
terial is more conducting than both equal-charge Y3Nb07 and Y3Nb07 though not as 
conducting as Zr2Y207. That equal-radius Y3Nb07 is more conducting than Y3Nb07 
itself is probably to be expected, what is more surprising is that it is also more con- 
ducting the equal-charge system which, as we have seen, has the kind of homogeneous 
distribution of mobilities seen in the highly conducting Zr2Y207 system. This is 
probably a consequence of using the Zr-0 short-range interaction parameters for the 
Y-0 and Nb-0 terms. As we mentioned above, these are softer than the ones for 
the Y-0 interaction but harder than the Nb-0 ones. However, since in this system 
there are three times more Y'^+ ions than Nb^+ ions, this means that, on average, the 
anion-cation short-range repulsion in this system is softer and this greatly enhances 
the anion mobility. A similar phenomenon was observed in the case of Ge02 where 
rescaling the short-range interaction terms by a factor of 1.6 enhances the oxide-ion 
diffusivity by several orders of magnitude |20j . 

As a final model system, we equalize both the cation radii and charges simul- 
taneously (i.e. we use the Zr-0 potential for all short-range interactions and set all 
cation charges equal to 3.5. As might be expected this gives a small distribution 
of tetrahedral volume sizes (figure [s] line with circles) and virtually no oxide ion 
trapping. This material is an incredibly good ionic conductor with a conductivity 
of ctme ~ 0.6 ri~^cm~^ at 1500 K, a value higher than the one found in the best 
yttria-stabilized zirconias (cr « 0.6 Q.-^CYn.-^ at 1667 K for 8% YSZ [21]). 



4. Vacancy ordering effects 

So far we have examined the relationship between disorder and mobility by examining 
the ionic positions and the effect of making particular changes in the interaction 
potential. In the Introduction, we motivated our study of this system by reference to 
the properties of the vacancies and the way they are ordered by their interactions with 
the ions. Their effect is, of course, implicit in the properties we have examined so far, 
but now we turn to an examination of the simulation trajectories by considering the 
vacancies explicitly. 

Vacancy interactions 

In [Appendix B| we explain how we identify the vacancies and how we can obtain the 
ion- vacancy and vacancy- vacancy radial distribution functions. Integration of these 
can be used to define coordination numbers. For example, integrating the Y'^+-vacancy 
rdf, gY~Y{r), from zero out to the position, Tc, of first minimum of the gY-o{f) gives 
the average number of vacancies in the first coordination shell of a Y'^^ ion. In figure [t] 
we show the cation- vacancy radial distribution functions for Y3Nb07 and equal-charge 
Y3Nb07 at 1500 K. From this figure the tendency of the vacancies to bind to 
ions is clear for Y3Nb07. This tendency is strongly reduced when the cation charges 
are equalized. 

In table |2] we report the number of vacancies around the different cations com- 
pared to the number which would be expected if the vacancies occupied sites randomly. 
It is clear from the table that vacancies bind very strongly to Nb^+ ions in rf-fluorite 
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Figure 7. Yttrium- Vacancy (solid line) and Niobium- Vacancy (dashed 
line) radial distribution functions in d-fluorite structured YgNbOj (top) 
and equal-ctiarge Y^NbOr (bottom). 



Table 2. Cation- Vacancy relative coordination numbers extracted from our MD 
simulations at 1500 K. 



Material 


ny-v 


nNb-V 


r^Zr-V 


RANDOM 


1 


1 


1 


YaNbOr 


0.81 


1.53 




ZraYaOr 


0.79 




1.21 


Equal-charge YsNbOy 


0.92 


1.23 




Partially ordered YaNbOr 


0.72 


1.85 




Pyrochlore Zr2Y207 


0.15 




1.85 



YsNbOy and less strongly to Zr'^^ in Zr2Y207. Considering the Coulombic effects, 
the vacancies (which have a relative charge of 2+) would be expected to bind more 
strongly to the cation with the lowest charge, i.e. Y^'^ in both materials. The fact 
that the contrary happens is therefore indicative of the importance of strain effects 
which are strong enough to outweigh the Coulombic interactions i.e. the vacancies 
prefer to bind to the smallest cations, despite their higher charges. This is confirmed 
by solid state NMR studies [22] on yttria-doped zirconias. 

We may now extend this analysis to examine the vacancy-vacancy ordering by 
integrating .gv-v between limits which bound the different coordination shells around 
an anion site in an ideal fluorite lattice. The first such shell gives the probability 
of finding a pair of vacancies as nearest-neighbours {i.e. along the (100) direction 
of the fluorite lattice) whereas the second shell relates to vacancy pairs along (110) 
and the third to (111). These probabilities may be compared with what would 
be expected if the vacancies are randomly distributed on the simple cubic anion 
lattice and we give the ratios in table [3] for simulations at 1500 K. As we note in 
the [Appendix B it is likely that the way that we identify vacancies leads to an 



overestimate of the number of vacancy nearest-neighbours. Nevertheless, it can be 
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Table 3. Relative population of the anion vacancy pairs expected from a random 
distribution of vacancies within a cubic fiuorite lattice, compared with those 
extracted from our MD simulations at 1500 K. 



Material 


(100> 


(110> 


(111) 


RANDOM 


0.231 


0.461 


0.308 


YsNbOr 


0.191 


0.418 


0.390 


ZraYaOr 


0.166 


0.391 


0.442 


Equal-charge YsNbO? 


0.136 


0.480 


0.384 


Partially ordered YsNbOr 


0.154 


0.304 


0.542 


Pyrochlore Zr2Y207 


0.073 


0.147 


0.780 



seen that there is a specific ordering tendency, with (111) vacancy pairs appearing 
with a higher frequency than would be expected for a random distribution, and 
nearest-neighbour pairs being lower. This order is in accord with Bogicevic's [3] 
examination of the energetics of vacancy ordering in lightly doped yttria-stabilised 
zirconia. The tendency we observe is, however, not as strong as might be expected 
from the calculations in this paper. This is, in part, due to the above-mentioned 
overestimate of the number of nearest-neighbours for algorithmic reasons, but is also 
because our calculations have been performed at high temperatures where the ordering 
is opposed by entropic effects. It is also interesting to notice that the vacancy ordering 
tendency decreases as temperature is raised and completely disappears at 2000 K. 
Finally, this tendency seems to be stronger in Zr2Y207 than in YsNbOr which might 
seem at first counterintuitive considering that Zr2Y207 is a better conductor than 
YsNbOy. This apparent contradiction will be explained in the next section. 

Effects oj vacancy interactions on the diffuse scattering 

The vacancy ordering influences the diffuse scattering observed in diffraction studies, 
and it is useful to compare the consequences of the degree of vacancy ordering we 
find in the simulations with those seen experimentally. Our simulations have repro- 
duced extremely well the structures derived from the analysis of the total (i.e. Bragg 
plus diffuse) powder diffraction [8] in these systems. However, several SAED studies 
[T5l [TB] have shown intense diffuse peaks for YaNbOy at particular positions in recip- 
rocal space. This suggests that the vacancy ordering effects are considerably stronger 
and of longer range than is indicated by the propensities we have shown in tables|2]and 
|3] In order to compare with these experimental studies we have calculated the diffuse 
scattering from the ionic configurations available from the simulations as described in 
[Appendix C| 

We decided to calculate the electron diffraction pattern of YaNbOr and Zr2Y207 
along the [110]/ zone axis (/ refers to the cubic fiuorite structure). This is found 
to be the most informative pattern in the experimental studies [HI [TBI IE!- In the 
experimental patterns of YsNbOy, a pair of intense and closely-spaced diffuse peaks 
is seen at Gp ± |(lll)*i where Gp is a reciprocal lattice vector of the fluorite-type 
lattice. In fact, the two peaks correspond to a ring of diffuse scattering centred at 
these positions which is projected down into the [110]/ plane. The Gp ± ^(lll)* po- 
sitions correspond to the location of Bragg peaks in the pyrochlore structure and the 
diffuse scattering seen in YaNbOy indicates a pattern of vacancy pairs along the [111]/ 
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direction which form strings along [110]/, as in the pyrochlore structure, but of finite 
range. [J] In Zr2Y207 the diffuse scattering is considerably weaker and the pattern of 
peaks is somewhat different. It appears as streaks along [001]/ in the studies by Irvine 
et al |16j and taken as indicative of C-phase ordering of vacancies (i.e. vacancy pairs 
along [110]/). Whittle et al [T7|, however, observe the pattern as 3 peaks centred at 
Gf± 5(111)* and interpret it as another type of pyrochlore-like modulated structure. 

In figure[8]we show the electron diffraction patterns calculated from our MD data 
(obtained from a simulation on 11,000 atoms, i.e. 10 x 10 x 10 unit cells). The pattern 
does not show the diffuse spots at Gf±5(111)* which have been observed in the exper- 
imental SAED patterns [151 IHl EZj • The calculated patterns show some diffuse signal 
at some of these positions but this appears weaker and less sharply peaked than in the 
experimental studies. This finding was confirmed in longer MD simulations on smaller 
cells. As a test, we made a similar comparison with the diffuse neutron scattering 
derived from single-crystal studies of (Zr02)i-x-(Y203)2; for 0.1 < a; < 0.25 and this 
was very successful, thus indicating that the disagreement with the experimental data 
is probably not a limitation of the way the diffuse scattering has been calculated. 

For YsNbOy we prepared a starting configuration for an MD run in which the 
cations were distributed randomly over the cation sublattice but in which the oxide 
ions were placed as if in a pyrochlore structure, to give the associated vacancy-ordered 
structure. In simulations at 1500 K, this initial configuration rapidly evolved to one 
in which the extent of vacancy ordering had relaxed to that seen in the previous d- 
fluorite simulations. It therefore does not appear that we can reproduce the evidence 
for vacancy ordering from the SAED studies with simulations based upon the d-fluorite 
structure with randomly disposed cations. 

In summary then, the rf-fluorite simulations reproduce the experimentally 
observed trends in disorder and conductivity and have enabled us to account for the 
way that these properties are affected by the cation size and charge. However, there 
is a suggestion that these simulations overestimate the conductivities, especially on 
the Nb-rich side of the composition range. Furthermore, although the simulations 
demonstrate vacancy ordering tendencies similar to those which have been found in 
more lightly doped zirconias these are not sufficiently strong to reproduce the SAED 
patterns. We are therefore led to reconsider one of the underlying assumptions of the 
work described so far, namely that the different cation species are randomly distributed 
throughout the lattice. 

5. Local cation ordering effects 

Our simulation studies of YsNbOy in the d-fluorite structure have shown lattice dis- 
tortions induced by the Nb^"*" ions, especially when two or more Nb^+ ions are found 

^ An ideal pyrochlore structure type, A2B2O7, is a superstructure of the fluorite (MX2) structure 
and is based upon a 2 x 2 x 2 unit cell with Fd3m symmetry (figures showing the differences between 
a fluorite and a pyrochlore structure can be found in refs ^ 1171 .^ The cation superstructure of 
pyrochlore is based upon ordering of A and B cations parallel to the (110) directions, separated 
by m , with respect to the origin. The A and B cations are, respectively, found at the 16c 
(eight-coordinated) and 16d (six-coordinated) sites whereas the anions are distributed between two 
tetrahedrally coordinated positions: 48f [0(1)] and 8a [0(2)]. There is another tetrahedral site 
potentially available for the anions, 8b, which is systematically vacant in fully ordered pyrochlores. 



Cation composition effects on oxide conductivity in the Zr2Y20-!-Y^NhO'! system. 16 
4 
3 
2 
' 1 





-1 1 

1 in [11m] 





-1 1 

1 in [lliti] 



Figure 8. Electron diffraction pattern obtained from the MD data for d- 
fluorite (top) and locally cation- ordered (bottom) Y^NbOj. 



as cation nearest-neighbours (the nearest neighbour Nb^^-Nb^'*' distance is 3.86 A as 
opposed to the ideal cation-cation distance in a fluorite structure a/y/{2) — 3.78A). 
The strength of this effect would suggest that the system will try to adapt by ordering 
the Nb^"*" ions so as to minimise the strain. Furthermore, the 0-Nb-O bond angle 
distributions extracted from the MD and RMC-generated configurations suggested a 
tendency towards octahedral local anion co-ordination of the Nb^"*" sites |S]. Moti- 
vated by these considerations and by our inability to reproduce the SAED patterns 
in YaNbOr we have been led to examine the consequences of local cation ordering in 
this system. Note that local (rather than long-ranged) ordering presents a significant 
problem for simulation studies because they necessarily involve periodic replication 
of a simulation cell which may be too small to allow for the decay of the structural 
correlations. Nevertheless, Miida et al [T5] estimate the correlation length of the order 
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responsible for the diffuse rings in the SAED pattern to be ~22 A which is the range 
of 4 fluorite unit cells (the majority of our d-fluorite simulations involved 4 unit cells 
along each cartesian direction, though the 11,000 atom MD cells used in the initial 
comparison with SAED involved 10). Furthermore, in the calculations to be described 
below, the cation ordering is imposed upon the system, using intuition and also the 
insight from the SAED studies to guess a structure which will minimise the lattice 
strain. In future work we will attempt to use a Monte Carlo cation-swapping algo- 
rithm to set up the system in a less biased way [21] . 

The natural stoichiometry for a pyrochlore-structured oxide is A2-B2O7 where B 
is a smaller cation which prefers octahedral coordination and A a larger one which 
prefers 8 oxide neighbours. We have seen above that there is a tendency for the 
smaller Nb5+ and ions to preferentially bind vacancies and hence reduce their 

coordination number relative to Y'^"'", and adoption of the pyrochlore arrangement 
would allow this tendency to be accommodated. Opposing this tendency may be 
long range strain effects, such as that noted above involving Nb^"*" in particular, and 
entropy, since the adoption of the pyrochlore structure involves ordering both the 
cation and anion sublattices with respect to rf-fluorite. Chao Jiang et al [18] estimate 
the configurational entropy as 

^Sideal = -4:kB[x In X + 

(1 - x) ln{l -x) + 2ylny + 2(1 - y) ln{l - y)] (3) 

where x is stoichiometric ratio of cations (| for Zr2Y207 and j for YaNbOr ) and y 
the fraction of vacancies to oxide ions in the unit cell | . This gives configurational 
entropies of 4.986 x 10"* eV/K and 4.534 x lO"'' eV/K for Zr2Y207 and Y3Nb07, 
respectively. Zr2Y207 has the appropriate stoichiometry to adopt the pyrochlore 
structure but the radius ratio of the cations is not sufficiently large to allow a stable, 
long-range ordered structure to form |18j - certainly not at the temperatures at which 
the solid-state synthesis is carried out. The situation must be very marginal however, 
as Hf2Y207 does form a pyrochlore [T7JIIH]) and the Hf'*+ ion is only slightly smaller 
than Zr''^"'" (note that Chao Jiang et al [T5] suggest that "covalent" effects may play 
a role alongside ion size in this comparison). Some local pyrochlore-like ordering in 
Zr2Y207 does therefore, seem plausible. The stoichiometry of Y3Nb07 is not that 
of a pyrochlore, but that does not rule out a pattern of local cation and vacancy 
ordering which is pyrochlore-like. Note that Irvine et al |Tnj consider cafion-ordering 
in Y3Nb07 likely and Miida [TS] has interpreted the SAED patterns by proposing 
modulated structures in which there is some local ordering of both vacancies and 
cations. The RMC-interpreted powder neutron studies [5] showed no evidence for 
cation ordering in either Zr2Y207 or Y3Nb07, but it must be noted that the neutron 
scattering lengths of Y'^^ , Zr^+ , and Nb^"*" are very similar and therefore no evidence 
of cation ordering can be obtained with this technique. 

5.1. Partially-ordered Y^NbO^ 

As mentioned above, Y3Nb07 does not have the correct stoichiometry to adopt a 
pyrochlore structure. However, it seems reasonable to assume that the smaller Nb^+ 
ions will prefer to occupy the octahedrally-coordinated B sites and we arranged the 
cations so that all the Nb^+ ions randomly occupy half of the B sites with the Y'^+ 
occupying the remaining B sites as well as all the A sites. As a consequence of the 
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random assignment of cations to the B sites the system does not exhibit "long" -range 
cation order within the simulation cell (though, because of the simulation boundary 
conditions, the system is periodic). The oxide ions were placed in the 0[1] and 0[2] 
positions. 

We ran a simulation on this material in the same way as explained above on a 
system with 11,000 atoms. The obtained electron diffraction spectrum is shown in 
figure^ It can be appreciated that this pattern shows peaks at Gp ± |(111)*, in 
reasonable agreement with the experimental data. The reason why we do not observe 
two diffuse spots but a single peak is that, because of the periodic boundary condi- 
tions, in reality we do have some long-range order in the cation positions. 

We computed the internal energies of the partially cation-ordered YsNbOy d- 
fluorite at low temperatures (T = 100 K). This was done by slowly cooling a well- 
equilibrated high temperature run. These runs contained 704 atoms (four fluorite unit 
cells in each direction). We find that the locally-ordered material has an energy which 
is lower by about /S.E « 56 kJ/mol = 0.58 eV/molecule than the d-fluorite mate- 
rial, confirming that some degree of cation ordering is indeed favoured on energetic 
grounds. From this we can calculate the order-disorder transition temperature, as in 
[IH], To-D ~ AE/AS,deai = 1160 K. This value goes up to Tq-d « 1460 K if we 
assume that only 2/3 of the vacancies are ordered, i.e. if we set y — 1/12 in equation 
|3] Considering that this material is usually synthesized at 1700-1800 K [8l[l5j|T6], it 
seems plausible that a certain degree of cation-ordering is observed depending upon 
how the sample is cooled back to ambient conditions. 

We find that the locally-ordered material is a much poorer ionic conductor than 
the disordered YaNbOy. We estimate that cation-ordered YsNbOr has a conductiv- 
ity of (7ne < 0.002 il^^cm^^ at 1500 K, approximately an order of magnitude lower 
than the d-fluorite simulations at the same temperature and in much better agreement 
with the experimental data. The reason for this much lower value of the conductiv- 
ity is that a partial cation ordering increases the Nb^^-vacancy association (see table 
|2| as well as the vacancy-vacancy pairing along the (111) direction (see table [3| and 
both factors hinder the vacancy mobility and therefore reduce the overall conductivity. 



5.2. Cation- ordering in Zr2Y2 J 

Cation ordering in Zr2Y207, on the other hand, seems to be a less important issue. 
We recall that the tendency for Zr''^ to preferentially bind vacancies, relative to Y'^"'", 
is considerably lower than for Nb^+ in YaNbOy, and the strain associated with the 
different local cation configurations is much less pronounced in this system. The ex- 
perimental SAED studies show slightly different diffuse scattering patterns (TB], [T7] . 
suggesting that lattice modulations in this material depend on the details of sample 
preparation. In all cases, the diffuse features are less intense than in YaNbOy. Whittle 
and co-workers [17] have estimated a correlation length of sa 11 A for the local order 

+ This value is found to depend weakly on the cooling rate and on the system history. 

* The error associated with this number is large because the corresponding mean squared 

displacement curve is very noisy due to the low mobility of the anions. 
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which is much smaller than in YaNbOy. 

Since Zr2Y207 has the correct stoichiometry for a pyrochlore, we have run 
simulations starting from a fully pyrochlore-ordered configuration (i.e. a long-ranged 
ordered structure). The internal energy of this system at 100 K is slightly lower 
[AE = 49.5 kJ/mol, a value in reasonable agreement with the trends observed in 
reference [18] ) than that of a rf-fluorite configuration obtained by cooling down a high 
temperature run and, considering that the ordered system has a lower entropy, it 
would seem that this structure is unlikely to become stable at higher temperatures 
(the estimated Tq-d is «1030 K). We also simulated the ordered system at 1500 
K to calculate the conductivity (note that there is no cation site exchange at this 
temperature, so that the system cannot relax to d-fluorite on the simulation timescale). 
We find the conductivity to be very small (at least two orders of magnitude lower 
than the experimental value) . An analysis of the vacancy- vacancy and cation- vacancy 
rdfs confirms that in pyrochlore structured Zr2Y207 there is a strong pairing of the 
vacancies along the (111) direction and that these bind very strongly to the Zr^+ 
cations (see table |2] and [s]) which explains the very small ionic conductivity shown by 
this sample. 

6. Summary and Conclusions 

We initially examined (i-fluorite simulations, with a disordered cation sublattice, which 
had been shown to reproduce very closely the neutron powder diffraction data [8] , and 
these already allow us to account for a substantial part of the experimentally observed 
trends in the conductivity and degree of disorder between YaNbOy and Zr2Y207. The 
difference between the two systems can be attributed to differences in cation size and 
to the different degree of strain imposed on the lattice by the highly charged cations - 
this creates oxide trapping sites in YsNbOy and gives an inhomogeneous character to 
the diffusive dynamics. Although both systems exhibit a tendency to order vacancies 
so that there is a likelihood of finding the vacancies close to the smaller cation, the 
vacancy ordering effects are not sufficiently strong to give rise to the diffuse scatter- 
ing features which have been observed in small-area electron diffraction studies. The 
fact that the observed trends are recovered from these d-fluorite simulations implies 
that they are primarily associated with the mean cation composition and traceable to 
the direct influence of ion size and charge, with any influence of intermediate-range 
vacancy structuring an additional feature. 

Because the d-fluorite simulations failed to reproduce the diffuse scattering seen 
in SAED studies we were led to consider the possible effects of partial cation ordering 
based upon the pattern suggested by the pyrochlore crystal structure. Our results 
in this context are only indicative as we have, as yet, only examined the properties 
of "by-hand" ordered systems. Nevertheless, in YaNbOr we were able to postulate 
a structure for which the diffuse scattering was in qualitative accord with the SAED 
diffuse scattering patterns and which was energetically favoured with respect to d- 
fluorite. In this structure the tendency of vacancies to bind to the smaller cation was 
considerably enhanced relative to d-fluorite and the conductivity was reduced by an 
order of magnitude, bringing it much closer to the experimental value. For Zr2Y207 
we were not able to postulate a partially ordered structure which would significantly 
improve upon the rf-fluorite calculations across the full range of observations (includ- 
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ing the SAED). In Zr2Y207 the calculated conductivity in the rf-fluorite structure is 
already close to the experimental value and the agreement with the powder diffraction 
is very good. The SAED studies indicate that the intermediate-range order is in any 
case weaker in this case and our own studies on d-fluorite show that the tendency for 
vacancies to associate with the smaller cation (the driving force for pyrochlore-like 
ordering) is considerably reduced relative to Y3Nb07. In this case a fully ordered py- 
rochlore is slightly favoured energetically over the d-fluorite structure but the entropy 
cost is high enough that the latter, more disordered structure, will be more stable 
above «1000 K. 

As we have indicated above, the finding that partial cation ordering has an 
influence on the material properties presents a challenge for MD simulation because 
at the temperatures at which we might hope to equilibrate the system with respect 
to cation positioning the ordering tendency has disappeared for entropic reasons. 
However, the cation ordering also presents a challenge for comparing experimental 
results. The materials are synthesised at high temperatures but then cooled (according 
to different protocols) in different studies: as the temperature is reduced the cation 
mobility will be drastically reduced so that the degree of cation ordering found in low 
temperature studies will depart from the thermodynamic equilibrium value to varying 
extents. 
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Appendix A. Simulation details 

We used dipole polarizable interaction potentials specified in reference f8^ (known as 
DIPPIM). These interaction models include full, formal charges on the ions and the 
Nb^"*", Zr'*''", Y'^+, O^^ ions have their full, in-crystal polarizabilities. The parameters 
for these potentials were obtained by "force-matching" them to ab-initio reference 
data as described in [23] . The full set of parameters is reported elsewhere [84 . Most of 
the simulations on Zro.5-o.52;Yo.5+o.25a;Nbo.25a;Oi.75 were done using a cubic simulation 
box with 256 cations and 448 oxygen ions (4x4x4 unit cells) though some calculations 
were done with larger cells. The time step used was (unless otherwise specified) 20 
a.u. — 4.84 X 10""* ps and all the runs were performed either in an NPT or NVT 
ensemble, with thermostats and barostats as described in reference [53]. Most of the 
high temperature runs were performed in an NVT ensemble, in which the cell volume 
was obtained from a previous run in an NPT ensemble with the external pressure set to 
zero. Coulombic and dispersion interactions were summed using Ewald summations 
while the short-range part of the potential was truncated to half the length of the 
simulation box, i.e. about 10 A . All the simulations were between 110 ps and 10 ns 
long. We used an MPI version of our code which we mainly ran on a small workstation 
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with 2 E5462 Xeon Processors. One such processor would yield about 500,000 steps 
{i.e. ~ 250 ps) in a day. 

Appendix B. Vacancy identification and calculation of the vacancy- cat ion 
coordination numbers 

We can identify the positions of vacancies for an instantaneous ionic configuration by 
finding which of the coordination tetrahedra are empty. Because the cations are not 
diffusing we can monitor the properties of each tetrahedron from the identities and 
instantaneous positions of the four cations which sit at its vertices. Such a tetrahe- 
dron is empty if no anion is within the volume bounded by the four planes defined by 
the positions of each set of three of the cations. Because the oxide ions may perform 
large- amplitude vibrations about their average sites at the temperatures of interest 
for dynamical studies we only assign a vacancy to a tetrahedral site if the tetrahedron 
has been empty for a minimum of two frames (i.e. 10 fs). This measure has already 
been used for the study of PbF2 [271 HI] and seems to identify vacancies reasonably 
faithfully even at high temperatures. For example, in simulations of Zr2Y207 at 1500 
K we find that the average number of vacancies is 86 whereas the stoichiometry would 
indicate 64. This number goes down to 67 at 300 K, in reasonable agreement with the 
stoichiometric value. Thermally excited Frenkel pairs may make a small contribution 
to the excess observed at high temperature but the most important factor is probably 
a shortcoming of our way of assigning vacancies. If a tetrahedral site contains a va- 
cancy, it will induce a substantial distortion of surrounding tetrahedra which lowers 
the reliability of the geometric criterion we use to decide if the site is filled or empty, 
especially at high temperatures where all the atoms are making large amplitude vi- 
brations. There is therefore a tendency to assign vacancies to sites which neighbour 
a site containing a true vacancy, and this is manifested in a larger than expected 
amplitude for the nearest-neighbour peak in the vacancy-vacancy radial distribution 
function discussed below. 

Once vacancy positions have been identified we can calculate radial distribution 
functions (rdfs) between ions and vacancies and also vacancy-vacancy radial 
distribution functions. Integration of the rdf can be used to define coordination 
numbers. For example integrating the Y^"'"-vacancy rdf, (?y-v from zero out to the 
position Tc of first minimum of the gy-o rdf gives the average number of vacancies in 
the first coordination shell of a Y^^ ion, ny-v 



where p is the density of vacancies. Analogously, integrating the vacancy-vacancy rdf, 
gv-y from zero out to the position Vc of first minimum of the (?v-v rdf gives the 
average number of vacancies in the first coordination shell of another vacancy. 

Appendix C. Diffuse scattering calculation 

We calculate the intensity of (total) scattering at a point q in reciprocal space from 




(B.l) 




(C.l) 



3 
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where the sums run over all ions in the sample, Vij is the distance between ions i and 
j and tti is the scattering amplitude of the species to which i belongs (for example, the 
neutron scattering length in the case of neutron scattering) . Because of the periodic 
boundary conditions, the accessible vectors q with a cubic simulation cell of side L 
are restricted to the set ^{m,n,p) where m, n, and p are integers, and this provides 
a limit to the resolution of the simulated pattern. 

We used the atomic position files obtained from the MD simulations on d-fluorite 
structures to calculate the diffuse scattering from equation |C.1[ The MD patterns 
were obtained from simulations on 11,000 atoms. These simulations were between 31 
and 100 ps long. We started them at 2000 K and the temperature was then lowered 
down to 300 K with a cooling rate lower than 50 K/ps. The and Nb^"*" atoms 
were placed at random on the cation lattice. We calculated /(q) for the six equivalent 
(110) planes and averaged in order to show the data in the figure. We input the values 
of /(q) on the grid of allowed q vectors into a contour program to generate the images 
shown in figure [8] 
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